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abstract 


Free surface problems are commonly met with 
in the fields of groundwater and hydraulic engineering, 
fef particular importance to Civil Engineer are the 
problems of flow into a well and flow over a spillway 
crest of a dam or a weir. These problems can be solved 
by various methods like graphical, model, analytical 
numerical methods. One of the numerical methods, 
finite element method has been used in the present study. 

I 

Finite element method is versatile and can easily 
incorporate nonhomogeneity , anisotropy and complex 
boundaries, ks such this method is preferred over the 
others, 

n 

The flow into a well in an unconfined aquifer is 
studied. A parametric study of the degree of anisotropy 
of the permeability has been carried out, A particular 
case of nonhomogeneity is also studied. In addition, the 
effect of partial penetration on the discharge and the 
free surface has been investigated. 

Experiments were conducted for the flow over a 
spillway using the Hele-Shaw model. It has been found 
that at low ratios of operating head to design head the 
model oannot predict the performance of the prototype as 
the viscous effects will be predominant * 



chapter 1 


INTRODUCTION 


The flow of water through a po^us media is xsf gr€!at 
importance in many fields of engineering, agriculture and 
geohydrology. Problems of unconfined flow are difficult to 
solve because of the exist ance of free surface or phreatic 
surface . 

Several attempts have been made to solve such 
problems and these may be classified as graphical, scale 
model, analogue, analytical and numerical methods. Graphical 
methods are generally used to solve steody state two-dimen- 
sional flow problems by construction of flow nets. The 
construction of flow net requires skill and is rather.te01ous.lt 
requires several trulls before one gets a flow net which 
satisfies all the boundary conditions. 

Scale model rand an&gue methods are an attempt to 
overcome the difficulties of a pureiy mathematical approach. 

Here the prototype flow through porous medium is replaced 
by a scale models Generally, however, it is difficult and 
time consuming to change from one configuration to another 
of proposed prototypes. Analogue models are subjected to 
similar differential equation and boundary conditions. 

Such analogues, however* are generally expensive to construct 
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and accurate measurements are difficult. In addition 
poor adaptability makes them difficult to change the 
model corresponding to a change in the design. 

Analytical solutions are normally possible only 

for those flow cases where boundary configuration and 
conditions are simple^ Consequently this le<i to certain 
simplified solutions of steady state two-dimensional 
seepage through homogeneous media. While the available 
solutions dj give reasonable results for certain 
situations > the idealized assumptions cannot be justified 
for many other cases. 

The numerical method of analysis is rapidly gaining 
in importance, especnlly with the advent of high speed 
digital computers. Till recently numerical approaches have 
generally used a finite difference solution of the governing 
differential equations. This method has been applied to 
a wide variety of seepage problems including confined and 
unconfined flow through porous media, 

A more recent method , namely , finite element method is 
,a numerical procedure which was originally developed for 
solving complex structural systems* Recently this method 
has been used successfully in the field of Soil Mechanics 
and Foundation Engine ering^ Rock Mechanics, Hydraulics and 
Water Resources Engineering, Nuclear Engineering, Mechanical 
Engineering, Biomedial Engineering, etc. 
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In the present study the finite element method is 
used to solve some free surface problems in groundwater 
flow hydraulics. The flow of v/ater into th e well in 
?n unconfined aquifer is stu^lied when the aquifer is 

t 

horn{j^£ieous and isotropic, h’^mogeneous and anisotropic 
for different percenb^igas nf drawdown. The case of non- 
homogeneity with isotropic and anisotropic conditions and 
partial penetration is also studied, 

Vascous flow analogue experiments were conducted 
to simulate the flow over a spillway. The flow over a 
spillway IS of highly contracting nature and such flow 
could bo assumed to be irrotational as long as there is 
no seperation. Such a flow problem is well suited to analysis 
by t he technique of potential theory and the results thus 
obtpine^ would be near to the aotu'al flow over the prototype 
structure . 

It IS well known that potential flow is well 
reproduced by menns of the flow of a viscous fluid between 
two parallel plates, sufficiently close to each other. 
Hele-Shaw model (or parallel plate model) was used for 
conducting the experiments. The results obtained have been 
compared with t he results obtained by various research 
workers. 



CHAPTER 2 


LITERATURE SURVEY 


The analysis of any two-dimensional ideal fluid, 
flow problem involving a free surface has attracted the 
attention of many research workers for a long time. Plow 
characteristics may include the velocity distribution or 
pressure distribution or both, over the entire flow field, 
and also the discharge coefficient, contraction coefficient 
and the free surface location. To this date exact solutions 
for these quantities have still not been found for many 
problems of practical interest. Most of the existing 
methods can solve only problems of simple geometric 
boundaries (1,2,3 end 4). 

Recent developments in finite element technique as 
applied to problems in various fields of engineering and 
the availability of high speed digital computer has attracted 
many research workers to solve some of the complex boundary 
Value problems. 

The first practical application of finite element 
method to fluid flow was to potential theory problems 
involving ground water seepage. The finite element method 
was first applied to seepage problems by Zienkiewicz (5 and 6), 
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Thereafter many groundwater 'problems have been solved 
by this technique, R.C. Taylor and G.B. Brown (7) have 

i 

considered the steady seepage flow vath a free surface and 
have given a trial and error method to obtain the correct 
free surface satisfying the boundary conditions. 

R,i3, Volker ( 8 ) has shown that the two commonly 
used head loss eouation for nonlinear flow in porous media 
will yield field equations which are amenable to solution 
using finite element approach, J.A, McCorquodale (9) has 
solved the non Darcy law using the finite element technique. 
He has shjwn that a C'jnverging successive over-relaxation 
solution for thr nonlinear elliptic partial differential 
equation for steady transitional free surface flow through 
porous media is possible when a variational approach is 
used. This paper deals with unsteady problems also 4 

Several workers have considered finite element 
solutions of the gener'^l fluid flow c^se e.g, the 
Navier-Stokes equation. J.T. Oden ( 10 ) has developed 
a finite element foimiulation dtsoribing a discrete model 
of compressible and incompressible fluids. He has shown 
that finite element '=inalogues describing the motion of 
compressible and incompressible fluid may be developed 
without resorting to variational principles by consider- 
ing energy balances. This paper is an outstanding example 
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of efforts to develop o, general finite element numerical 
solution to t he gener?al flow equations, 

P.F, Lemeux and T 4 E, Unny (11) have formulated 
a variational principle for the flow of a viscous inccmp*- 
ressibl*-^ fluid, which includes the convective terns and 
covers both timo-independent and time -dependent phenomena. 

J.A.* McGorquodalcand M,K. G-iratalla (12) have 
modified the two-dimensional conformal mapuing solution 
for the drag characteristics for a supercritical flow 
over a continuous sill to predict the drag on a dentated 
sill. They observed that t he finite element method hod 
the disadvantage of requiring a good approximate starting 
solution in order to avoid stability problems. The finite 
element solutions that were obtained were in good agreement 
with the experimental and theoretical drag coefficients. 

Liam Finn (13) has presented a s'^lutJon of two- 
dimensional steady state seepage problem through dams 
using a networh of linear triangular elements. A. free 
surface was initially assumed and subsequently adjusted 
untill free boundary c ond it ions were satisfied. The results 
compare well v/itb t he results obtained by other research 
workers. 

P.W. France et al (14) have e xtended the finite 
element discretization technique to transient three-dimensional 
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seepage problems bounded by a phreatic surface. A comparison 
of the numerical and analogue fSee surface configurations 
shows good agreenent , 

M, Ikegawa and K. Washlzu ( 15 ) have analyzed the 
flow over a spillway crest using variatunal principle 
an assumed free surface to begin with, and subsequent 
Iterative solution by the use of finite element technique. 



CHAPTER 3 


FORMULATION OF THE PROBLEM BY FINITE ELBMEHT METHOD 

3.1 INTRODUCTION TO FINITE FLEi^aNT METHOD; 

Recent developments in the finite element methods 
of solution can be interpreted interms of variational 
procedures (as the minimization of the total potential 
energy functional) leads naturally to its extension 
to other bounU ary v alue problems. In many such 
problems two equivalent, alternative, formulations 
exists. In the first a partial differential equation 
is v/ritten and its direct solution attempted. In the 
second, problem is to find a function minimizing a 
certain specified functional over the whole field 
involved. The two appro^^ches are mathematically equal, 
an exact solution of one being the solution of the other. 

The method of finite element dhvolves the' follow- 
ing steps in solving a problem. The domain of interest 
18 subdivided into a number of discrete elements, called 
finite elements. These are in general of polygonal 
shape. The unknown function over an element is assumed 
to vary in the form of a polynomial of the loc^l 
coordinates. In choosing a particular polynomial 
approximation certain convergence criterion has to be 
satisfied (discussed later in this chapter). The 
functional minimization is carried out for each element 
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over the whole region jby treating the nodal values as 
paremeters* This results in a senes of siiault aneouo 
equations* The solution of these equations gives an 
approximate solution of the problem. 

It IS of interest to remprk thnt while the well 
knovm finite difference method represents a direct 
approximation appronch to the differential equation 
type of formulation, the finite eloirent method is an 
approximation applied to the variational formulation. 

The freedom of choice of the shape of the element, and 
the case of treatment of nonhonogeneous situations make 
the finite element Tvig-thod much more versatile. 

Let the physical (or purely mathematical) formu- 
lation of the problem require the minimization of a 
functional E over a certain domain. E is defined as an 
integral over the domain T and part of its boundary S 
on which the unknown function, |h^ , or its derivatives 

appear, that is^ 

B = / f ( {h> , ijE ^ 

+ / g' (|h} , {h] , ) dS (5'.1) 

s 
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Let the region be divided into smaller parts, 
subregions triangulp.r in shape which we shall call 
elements, and let the function which we are trying to 
determine be described in each element as 


■jn} = [u] |h]® (3,2) 

In this ^ Hj may list nodal values of the 
function associated with the element or simply some 
parameters characteristic of it, Curly brackets are 
put round the unknown function to show that it may be 
a vector, N is shape function which is function 
of the coordinates only. 


To miniFiize the functional E with r espect to the 
total number of parameters {hj associated with t he 
whole domain we must write a system of equations 



aE 

m. 

j M 
‘ an- 


= 0 


. aE . 

V5H,/ 


( 3 . 3 ) 


If it is true to say .that the total functional 
IS equal to the sun of the contributions of each elemenb, 
that IS, 
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E = E E'^ 

then a typicel equation becomes 

as ^ ^ a^ 
aH^ 

With summation over all elements. The rule for 
bhe assembly of a general minimizing set of equations 
IS thns available. 

In the special case where E is a quadratic funct- 
ions. of and its derivotives we shoU find that we 

can always write the element derivatives as 

"wV - (3.6) 

- e e 

In which and are constants 

representing the element permeability matrix and nodal 
vecror of applied flow, respectively, Now the minimizing 
set of equations (3.3) can be simply written as 


(3.4) 


(3.5) 


- W Ih; Hi') - 0 

In which 

[k3 = E [kJ ® 

= £ -ii'} ^ 


(3.7) 


(3.8) 


(3.9) 
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3 . 2 CQNVj]-HGE])rCE CRITERIA 

The function approximation given by equation 
(3il) must obey certain completeness rules to generate' 
the conyergence of the results the subdivision into 
even smaller elements is attempted, 

First, as the elej ent size decreases, the functions 
f and of the integral (3.1) must tend to be single 
valued and well behave 1 in physical problems. Thus it ^is 
necessary that the following criterion should be satisfied. 

Criterion 1 ; 

The element shape functions [Nj must be such 
that With a suitable choice <Df [hSt any constant 
values of or its derivatives present in the functional 

E shuuldbe able to bo represented in the limit a s element 
size decreases to zero. 

Second, the validity of summation implie^i in 
equation (3.4) musm be preserved, and unless special 
integrals are aided for the inter-elenent boundaries 
we must make sure that the terms like f and remain 
finite there. This is achieved if highest derivatives 
of H occurring in these expressions are finite which 
leads to criterjion 2; The element shape funcbions [nI 
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have to be so ^hoosen that at elemeiit interfaces 

an(3 its 'derivatives aare of one order, lea-B than that 

occurrixoi in expressions f which define the 

u* 

functional, are contirona*. . 

3*3 ftJ)VAidTA&BS 4ND DISADVlHTA&SS • 

Like all nuiTierical approximations, the finite 
elemeit method is based on bhe concept of discretization. 

Not only (3oe3 the idealization portray the body as 
continous, but it also requires no soperate interpolation 
process to extend the approximate solution to every point 
within the continuum* Despite the fact that a solution 
IS obtained at a finite number of discrete nodal points, 
the formulation of the field variable models inherently 
provides c solution at all other locations in the body. 

In contrast to other variational and residual approaches, 
the finite element method does not require trial solutions 
which must all apply to the entire multidimensional continuum. 

Some of the most important advantages of the 
finite element method derive from the techniques of 
introducing boundary conditions. This is another area 
in which the method differs from other variational or 
residual approaclies. The introduction of boundary con(iitl'fl>ns 
into the assembled equations is relatively an process. 
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It IS simplified in that only the geometric boundary 
conditions need be specified in a variational approach 

because the natural conditions are implicitly satisfied. 

a 

The finite element method not only accomdates 
complex geometry and boundory c inditions, but it also has 
proven successful in representing various types of compli- 
cated material properties, that arc difficult to incorporate 
into other numerical methods* (Por example, formulations 
in seepage problems have beon devised for anisotropic, 
nonhomogeneity conditions). One of the most difficult 
problems encountered in applying numerical procedures of 
engincjering analysis is the represent a bion of nonhomogeneous 
continue. Nevertheless, the finite element mefhe^d readily 
accounts for nonhomo gene it y by the simple tactic of 
assigning different properties to different elements, If 
a refined represent^^tioii op the variation of the material 
characteristics is desired, it is even possible to vary 
the properties within the element according to a 
preselected polynomial pattern. The systematic generality 
of the finite element procedure makes it a powerful and 
versatile tool for a wide range of problems, ' 

' T 

Even the most efficient finite element computer 
program requires a relatively large amount of computer 
memory and time. Hence, use of the method is limited to 
those who have access to relatively l^rge, high-speed 
computers* 
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3 . 4 FORMUL/ITIQN Off THi) PRQBLt^M ; 


Th3 'Darcy's law of Seepage' establishes a linear 
relationship between the seepage velocity and hydraulic 
gradient, This law, which is a simple consequence of 
viscous flow neglecting inertia effects, can obviousl;^ 
be generalized to a tv;o dimensional situation. 

If the velocities in the direction of two 
orthogonal axes x and y are designated by a vector, V or 



(3.10) 


and if the head gradient is defined similarly by its two 
components, 



(3.11) 


Then the most general line^'^r relationship that can exist 
between the two quantities is of the form 


{vi = - [k] {grad H| 


( 3 . 12 ) 


In which is a 2 by 2 matrix of the coefficient of 

permeabilitv , le. 


[k] = 


K 


XX 


0 


K 


yy 


0 
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If the dirtiction of the axes is changed to x' and 
y', then the velocity vector in the direction of new 
axes can be found as 


{v'j = i l] = [lI iv) (3.13) 

In which '.ij is the transform. s'! ion matrix of directional 
cosines. Similarly the new vector of head grac^ient 


-.'grad H’i = 

1 j 


aH’l 


9x 



= - [ijj {grad H? (5.14) 


Combining equatiois 3.12, 3.13 and 3.14 yields# 

(v'j = - [k*] {graa H’} (3.15) 

In which the nev/ permeability is 

[ic] = [l] [k] llT’' (3.16) 


V/ith this tyi'e of trensf ornation, and if K is symmetrical, 
it IS always possible to two orthogonal directions 

for which K’ reduces to a diagonal matrix giving 











(3.17) 
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These directions are known as the principal axes 
of the porous material and the particularly simple 
relationship, eq.3.15,in the direction of these axes 
are obviously worth noting* In a stratified material 
such as rock or soil the behaviour in a direction perpendi- 
cular to the strata is symmetrical and this axis coincides 
with one of the principal direction. 


3 . 5 SOLUTION OF TWO-DIjVTEdSIQNAL SEEPA&E PROBLEMS 

If the X and y directions (locally) coincides with 
the principal axes of the material, then the continuity 
rel'^'tion 


6u ^ 

dx ay 


0 


(3.18) 


Changes by substitution of eq.3.17 into the 
governing differential equation 


Avhich IS valid for both homogeneous and nonhomogeneous 
situations. 

The solution of eq, 3.19 subject to specified 
boundary conditions specifies the problem in a unique 
manner. However, an alternat^ formulation is possible 
with the aid of the calculus of variation. The well known 
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Euler theorem states th'^t if the integr^'l 


E(u) = fff f( x,y,H II , ||-) dx dy 


(3.19a) 


IS to be minimncd over a bounded region V, then the 
necessary and sufficient condition for this minimum to 
be reached is that the unknown function H (x,y) should 
satisfy the following differential equation, 


a j af I , a 

^ (^FTWVT/ ^ 


J 01 . 

1 a'(aH/ay) J 


af 

an 


= 0 


(3.19b) 


With in the same region, provided H satisfies the sarae 
bound=^ry condition in both the enses. 


It can be observed that the equivalent formulation 
to that of eq, 5i19/tho requirement that the volume 
integral given below and taken over the whole region 
shoulc? be minimized. 


E 


Iff 

V 




an. 2 


+ K. 



dx dy 

( 3 . 20 ) 


subject to H obeying the same boundary condition. 

This IS valid whether K and K are constant or 

,y .y 

variable with the x and y coordinates. Hov/ever it is 
possible to give it a physical interpretation. The expression 
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V. 


IS 

0X 


+ V. 


y 


dH ^ 


- K 


- K 


(5.21) 


represents the rate at which energy iS being dissipated 
in a unit volume. The integral of eq, 5.21 is proportional 
to the rate of energy dissipation over the entire region 
and the fact that a pressure system which coincides 
with a minimuni of this gives the true solution is 
consistent with the v/ell-lmown principle of 'minimum work’. 
The problem studied is shown in Fig, 1. 

Let the region in question be divided into 
arbitrnry trlangul^^r areas or finite laments, as shown m 
Pig. 2 andj3, IT the unlcnown values of the function H, at 
the nodes of each element define the function throughout 
the entire region completely and uniquely, then different- 
ic'bing the functional E with respect to each of these 
nodal values and equating each of the differential to zero 

ijt 


X 



Fig. 2 
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will result in a series of airnultaneous equqtions< 

From these the final solution can be obtained by 
feeding in the appropriate boundary values. 

Considering a typical triangle i»3,m with a 
locally defined set of coordinates of the nodes^the 
simplest definition of the function H v/ith the element 
will be obtained by talcing a linear function 

H = g + bx + cy (3*22) 

where a, b and c are constants. 

Evaluating the three constants interms of the 
coordinates of the nodes an- the corresponding nodal 
values of the function. Perfoiming the necessary algebra 
(in matrix language) yields 

(3.23) 

(3.24) 

function 
Th^ coefficients 


H = |^(a^ + b^x + c^y), (c^^ + b^x f c^v), 

or simply 
H = 

in whioh{H®} stands for the values of the 
characteristic of the element considered. 


r I 


H 


H, 


H 


are defined as 
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a 


1 


b 

1 


2A 

- y,i 
2 A 


c 


1 



(3*25) 


With :>theTs followi ig a cyclic, anticlockwise, order 
in 1 , j, and In eq. 3.25 A is "the area of the 
trian'^le i, ra. 


To obrqin the ninimmr of the function E, it is 
necessary only to c^erive the necessary conditions for 
each cloiaent s^psirately. The influence of other elements 
on the (general mini^ium equations will follow in identical 
pattern, and to obtain the i^^neral equation we add up the 
contrj buto on of each element in the whole region. 


It 13 ^ importance to note here that an abrupt 
varnatijn of properties between elements defined by the 
vr\lue of If in equation 3*19 and 3.20 docs not introduce 
any difficulty as in eqn, 3.20 no infinite terns amse. 

Thus termiiiq E^, the contributun of the element 
i,j,n yields 



If 


K 


aH 3 


XX 3x 




(i) 


+ K 


yy 


9H 


5 


(3.26) 


dx dy 


22 


I 


With the integration limited to the area of the 
triangle i, j, m. This on substitution of eqn« 3.24 
becomes 


i: =//K 


XX 


dN^ 

5x" ’ 'Cx~ 


{h j (•^) dx dy 

OX 


r 


//K 


yy 


1 j mi 

dy * * Oy 


|h I ( ^ ) dx 


dy 


(3.27) 


Recognizing that the element 'contributes only to 
the differentials with respect tu the values at its three 
nodal points, assuming the permeability coefficients 
K and K to be constant with in the element, substitution 

Jr y 

of eq. 3.25 yields simply 



In which the [sj matrix hf^s coefficients of the type 




(3.29) 
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Tu assemble a typical equation for a differential 
of E with respect to any nodal value for the whole region, 
the contributions of elements adjacent to the node 
only are nonzero and the cssembly will be clearly according 
to the pattern, 




(3.30) 


or 


I T. S 


n 


in 


= 0 


(3.31 ) 


In the eqn. 3.31 is Icnown and hence the 
potential H at each node can be calculated. Knowing 
H, the velocity and hence the discharge bhrough each 
element can be calculated by eqn. 3.12, 
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3.6 BOUND^IRY CQttoITIONS 

In the solution of seepage problems three types 
of boundory conditions m^y he specified, i) surfaces of 
Known pressure (t^ead), li) surface of Known flow and 
111 ) iMiscible interlaces between two fluids. The 
v^rintion of pressure, 6p , on o surface is if the 

nodal pressures are set eaual to the boundary pressure 
octinr” at the no He, On the other h^nd, known conditions 
of flow "*t surfaces -^re handled by conputing the equivalent 
nodal flow in each element and then 

The interface between two inmisible fluids is 
'"iescrib-''d ’ry the two continuity equation as 

p'' = (3.31) 

'xlong Cj an-i <ln “ '^n (3.32) 

In which Cj is the interface, subscript 1 and 2 
refers to fluid 1 and 2 and n is a normal to the inter- 
face surface and p and q are nodal pressures and flow, 
respectively. In finite element f ormulation eqn. i3.31 
is satisfied only at the nodal points. Along the interface 
eqn. 3, 32 is dealt intorms of nodal flow contributions which 
account for the tributary regions, about nodal points 
alone: the interface. While this statement of the interface 
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GOnliti'^ns IS forw'^rl, tho solutj-on of a 

specific pron]t^ is complicate! by th'-* necassity of 
knvjvan^ ttic position of Cj. In the solution of these 
problems, it is nccoss'^ry to locotc while s-itisfyins 
eqns. 3.51 and 3.32, 

To lLmon>strate , c oisidur 3 . free surface. The ?\ir 
flow is ne^^iected oncl eq. 3.31 and 3-32 can he written as 

P - 0 on Cj ^ (3.33) 

and q = 0 ( 3 . 34 ) 

In which aero pressure corresponds to atmospheric 
(reference) pressure. 

In aatisfvinq bound ^^ry conditi ms in finite element 
method, it is possible to specify either uressi^re or flow 
rates, bub not both siraultaneously . For a free a urface 
these are three conditions to determine or satisfy na 
follows : 

(1) Determine Cj- (2) n == 0 and (3) = 0. 

Because it is not in .general, passible to satisfy 
conditions ( 2 ) and ( 3 ) simultaneously, it is necessary to 
assume a location of C-j-, set either condition 2 or 3 and 
subsequently adjust Cj unbil^ both conditions 2 and 3 are 
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sntisfiud. In this work the order of opero.tion. of this 
essentnlly variationol pppronch is the initial 

locetinn of Cj (b) = 0 on Cj and (c ) investi^^ation 

of p with subsequent relocation of Cj an(^ repetition until^ 
p 0, This oethoi was selected because of its sinplicity. 
The anibi -rui in th*- st ocification of these conditions 
must e^ist, where th^ free surface, Cj, intersects the 
/geometric boundary of the l^ody, because of nodal flow 
at this point can not be zero. This difficulty is 
miniraa zed by reducin/r the mesh size adjacent to that 
point. 

First a possible free surface is assuned , 

Wo flow IS allowed across the surface, iiid the finite 

element solution is carried out. The resulting pressures 

indicate tho direction an^i m>£mtude of the ch"^n»^e of frea 

surface that is necessary to satisfy equation 5 and 4* 

2 

Such a chans^e in the free surface is made to define Cj- 
by novin:^ each node along prescribed directions (v/here 
naterial b^nd exists, these directions coincide with 
their interfaces). li new finite element solution is 
obtainel and the reSultin^r n'>n-homogeneitiea of pressure 
indicate the steps necessary to locate C£, The procedure 
IB repeated untilffi the free surface pressures all approach 


zero . 
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V ISCOUS TLOV/ MOTEL STJDIES 

Oroundv/ater flow problems can be studied by 
Various means, Analytical methods c an be used where the 
boundary conditions nre simple, ’j/here the analytical 
metbiods become too difficult, recourse is made to numerical 
methods '^nd/or model studios. 

Ground water models can be grouped into four 
general types. S'Uid, electrical, viscous fluid and 
membrane. In the present studv it is proposed to study 
one of the above models, n^unely , viscous flow model or 
parallel plate moOel, This is clso called Hele-Shaw 
model after Helc-Shaw. 

4 , 1 DESCRIPTION Of THE MOTEIi; 

For the steady-state flow of viscous, incompressiblo 
fluids (at small Reynolds number), the Navier-Stokes 
equation of motion, the most general equation governing 
fluid flow, reduces in form to generalized statements of 
Darcy's lav/. Recongnising this Hele-Shaw in 1897 devised 
an appar^ntus whereby two dimensional graundweter flow 
could b3 investigated experimentally for structures with 
complex boundaries (3t 16). A viscous fluid is allowed 
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to fljv^ between two cloaely sp=iced parallel plates* Glyserine 
is used in the present study. 

The HDdel consists of tv/o perspex plates 1,3 cm 
thick set on => base plate v/ith the help of angles. 

Spacers are provided between the perspex plates in the 
form of washers. The plates are 60 cm long and 30 cm high. 

At the end of the plates ther'e is a constant head 
reservoir 8 cm. dia and 60 cm. high connected to a over- 
head reservoir 11 cm, dia. and 30 cm. high. The rate of 
inflow into the model can be adjusted, A graduated 
jar is kept at t he other end to measure the outflow (I’ig* 5). The 
shape of the structure to be studied is kept between 
the plates (17). 

4*2 VISCOUS FLOW MODEL THEORY : 

The Navi er- Stokes equation of motion, the mast 
general equations governing fluid flow, maj/’ be expressed 
as, 
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where 


X, y, z = spatial coordinates 

u, V, w = velocity comp)oneits in x, y, z directions 
respectively. 


X, Y, Z 



= body forces in x, y, z rUrections 
respectively . 




mass density of the fluid 
-specific vveight of the liquid 
g-acceleration due to gravity 




kinematic viscosity 
dynamic viscosity 


D 



+ V 


8 

W 




w 


a 

8z 


2 m 

If the distance/between the plates is small, 
the flow becomes two-dimensional (w=0). In rectangular 
coordinates the x-axis will be choosen as horizontal 
and midway between the plates, the y-axis vertical and 
the z-axis perpendiculor to the plates (Fig. 4). Because 
the only body forces are due to gravity, assuLiing a 

gravity potential of/'- = gH, where H is the elevation 
head (H = h - p/f g), we have 

X = - , Y = - and Z = - 

8x 6y 


3(gH) 

9z 
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For incompressible fluids, the third terms in equotion 1 
4*1, ywhich represents volume chpnge, vanishes. Since the 
velocity et the pl^ites must be zero (Fig. 4 b) the change 
in velocity components (a end v) with respect to z will 
be much greater than the velocity changes with respect to 
the X and y directions. Hence, in comp am sis'"' to their 
respective derivatives in t he z direction, we may assume 

, 0u , 9v and 4^ , and their second derivatives 

^ 

are of negligible order, How assuming steady-state 
conditions and making the stated assumptions, equation 4.1 
reduces bo 


^ ( *? gH + p ) _ u 

9x 




Vv 


9 ( ?gH + p) 

3y 




( 9 gH + p ) ^ Q 

9z 


(4.2) 
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\z (Pig. 4b) 
Fig. 4 Sketch showing parallel plates Model 
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The third of equation 4.2 shows thcit the total head 
at any point within the flow domain will depend upon 
X and y coordinates only. Hence the first two equations 
be integrated with respect to z to yield. 


Z 


c'x ' f]z 


Z 


9h _ .y vlv 

?y r az 


where the constants of integration are zero (from symmetry). 


= 0 a,t z = 0. Integrating once more 

with respect to z and noting that u = v= 0atz = + m. 
V/e obtain ■for the velocity components. 


u 


V 


+ ~ ° — ( f gH + p) 

2 f'* Px 


(4.3) 


Which demonstrate a parabolic velocity distribution 
(Pig. 4b). The maximum velocity at the centre of the 
channel will be (at z = 0) 


u 


V 


Sif ^ gH + p) 

2 /'A 9x 




( fgH + p) 


(4.4) 


and the mean velocity will be 2/3 the maximu^ veJ.cjC'^tjf/.y) Mrijfi 

CENTRAL lISIx'iRY 
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jji 


J1 


( + p) 


0:c 



0 

( fgH + 

p) 

(4.5) 

Multiplying^ and dividing by 


in the RHS ^nd setting 

H + -E- 
) g 

= h we 

got , 



'V. = - 

m 1 g 

at. 



3 ft 




V = - 


Qh 


(4.6) 

m 

3/'< 

0y 




‘‘/e see the nn‘=tlog 7 with Darcy's law. Darcy's law can 
be ejcprcssed as 


V = - KP II 

? g Sh 
" 3f‘ es 

of 

where KPis coeff icient/pernieabili by and 


( 4.7) 


s stands for coordinate direction along the stream- 
line , g, A/ have the same meaning as before. 


4.3 SCALE RATIO: 

Generally physical models are constructed using 
suitable scale ratios in true and distorted models to 
study the behaviour of prototypes. In the present analysis 
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Reynolcis number vail £PVern flow. The followins scale 
ratios pre ieriv-d asim'^ Reynolds law. 


(I?) 

^ p ^ m 


7 P 


Where rn and p refers to model and prototype respect' 
ively. V is the velocity, Ltis the linear dimension and 
th^ kinematic viscosity of the liquid. The velocity 
scale ratio, can bo derived as 


V. 


V, 


m 


V 






T iFie scale rat^o is not considered here since 
only steady st ^t e pro bleni i s analysed , 
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DISCUSSION OF RESULTS 


Cjmput^-cions fur the follJwing cases have been 
carried out using the Finite Element Method, 

^ • FuIIy Penetrating Well 

(a) Homogeneous and IsotropJC 

(b) Homogeneous and (\nisotrop;yf 

(c) Non-homogeneous and InisotropyC' 

2 • Partially Penetrating V/ell 

Homogeneous and Isotropic 

1 . Fully Pentrating Well 

A finite element mesh consisting of 126 nodes and 104- 
elements is used for the computation (Fig, 3). The free 
surface coordinates satisfying the boundary conditions and 
the element flow is obtained. The results obtained are 
shown in figures 6 to 11, 

( a ) Homogeneous and Isotroptd’ 

The free surface is obtained for different 
percentages of drawdown, namely 25, 37*5 , 50, 625 and 75* 

The representative free surface with respect to distance 
from the well is shown in Fig, 6 a and b. From the figure 
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It can be seen that the free surface obtaine-l by the finite 
element method is above that of an^^lytical solution near 
the well and the percentage difference between these t wo free 
surfaces v?ries from 0 to 5 per cent gradually increasing with 
the decrease of drawdovm. Beyond 20' fro^i the v/ell the values 
more ior less couicids, the naximura difference that is 
observed is only of the order 1*5 per cent. 

b , Honogeneous anl \nisotrop:^0 

In fig. 7. a 13 shown the variation of 
relative discharge with respect to anisotropy of 

the aquifer expressed as a ratio of permeabilities (K.^/K. 2 ) 
for different drawdown. It is observed that the variation is 
linear for all percentages of drawdown. The relative discharge 
decreases as the percentage of drawdown is increased for a 
particular va.lue of K^/K .2 and for the sane percentage of 
drav/d:>wn, the discharge increases with the increase in 
anisotropic ratio. 

The drawdown with respect to x, the distance from the 
v/ell is shown in figures 8, 9 and 10 for different values 
of anisotropic ratio (K^/Kg). It is seen that the point of 
inter section of the free surface at the well increases 
as the anisotropic rabio is increased i.e^the discharge 
increases as the anisotropic ratio is increased. 


36 


c ) Nonho mo^eneous end AnisotropTd 

\ particular case of nonhomogeneity of the lower 
layer having twice the pemo<ability of the upper layer 
IS considered for the purpose of illustration* In 
each layer different '^e^.rees of anisotropy are considered, 
Pig, 7.b shov;s the variation of relative discharge 

1. e* with respect to anisotropic ratio. 

It IS observed that the discharge in the nonhon ogeneous 
aquifer is relatively large when compared with homogeneous 
media for a given anisotropic ratio. In fig, 11--a the 
values of drawdown are plotted against the distance 
from the v/ell for different permeabilities. Prom this 
it is egain verified th'^t drawdown at the well increases 
with the increase in the permeability, 

2 , Partially Penetrating V/ell; 

Discharges were computed for different percentages 
of well penetration. A graph of Qj>/Qp with respoct 
to the well penetr?tion is shown in Pig, H-^b^ The 
values obtained are compared with the one given by 
Muscat (1), It is observed that the graph is similar 
to that of Muskat, but slightly ebove it. This can be 
explained by the fact that the discharge varies with the 
R/rw ratio. Muskat has taken this as 2000 where as 
in the present study it is 40, Bereli (20) has verified 
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TJ 

this in his paper* He has taken ~ = 112 for which 

J;W 

discharge computed vanes appreciably from that of 
Mu skat . 

F low Over Spillway: 

The discharge coefficient, 0^ as a function of 

hg/hp IS shown in Fig. 12. It is observed that the 

tig 

maximum discharge coefficient occurs st ^ = 0*67 and 

0.95 for P = 6'' and 4*' respectively. After it reaches 
maximum reduces and will be constant after a certain 
(hg/hp). The result do not conform with that of 
prototype. This may be due to the fact that viscous 
forces are predominant in the present case where as the 
gravitational effects are predominant in the prototype. 

This fact has been reported by Cassidy (21). As the height 
of the spillway) P, is decreased it is observed that t he 
maximum discharge coefficient will occur at a higher value 
of hs/hp. 

The present work is a parametric study with the 
use of the finite element method as applied to ground 
water problems. The results obtained are encouraging. The 
same computer program could be used to get results for 
gome of the difficult problems. 
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